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! Abstract 

>: 

{Sj The rate of electron-positron pair production in linearly polarized counter-propagating lasers is 

evaluated from a recently discovered solution of the time-dependent Dirac equation. The latter is 

^: 

■ solved in momentum space where it is formally equivalent to the Schrodinger equation describing 

o 

£Sj , a strongly driven two-level system. The solution is found from a simple transformation of the 

Dirac equation and is given in compact form in terms of the doubly-confluent Heun's function. 

• rH ■ 

^ ■ By using the analogy with the two-level system, it is shown that for high-intensity lasers, pair 

production occurs through periodic non-adiabatic transitions when the adiabatic energy gap is 
minimal. These transitions give rise to an intricate interference pattern in the pair spectrum, 
reminiscent of the Landau-Zener-Stiickelberg phenomenon in molecular physics: the accumulated 
phase result in constructive or destructive interference. The adiabatic-impulse model is used to 
study this phenomenon and shows an excellent agreement with the exact result. 
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I. INTRODUCTION 



The production of electron-positron pairs from classical external fields has a long his- 
tory, starting with the seminal work of Schwinger , where antimatter production from 
a constant electric field was considered. On the theoretical side, Schwinger's mechanism 
is relatively well understood and is usually interpreted as the decay of the vacuum into a 
particle-antiparticle pairs (in the Dirac sea picture, this is seen a a tunneling from the neg- 
ative to the positive energy sea). However, an experimental validation of this phenomenon 
is still out of reach: the intensity of a laser electric field required to produce an observable 
amount of pairs is on the order of 10 29 W/cm 2 [3j, which is still unattainable experimentally. 
This occurs because Schwinger's result states that the probability to produce a pair (per 



unit volume and time) results from a tunnelling process and is given by jl 



2] 



P s ~ e N s 



(1) 



where m is the electron mass, c the speed of light, |e| the absolute value of the electron charge 
and E the electric field strength. Thus, an appreciable amount of pairs can be produced 
only if E ~ E s = ~ 10 18 V/m, which is much larger than Ei s = 5 x 10 11 V/m, the 
typical electric field in the ground state of an atom (Is orbital). 

In the last few decades, laser technologies have made a giant leap forward such that elec- 
tric field of unprecedented intensity level can be attained (on the order of 10 22 W/cm 2 and 
higher 4j). At these intensities, relativistic effects start to be important and thus, the spon- 
taneous creation of electron-positron pairs from laser fields becomes more plausible. This 
has triggered many theoretical studies recently where pair production from generalizations 
or variations of Schwinger's process were considered. Thus, different field configurations to 
produce pairs have been studied such as counterpropagating lasers [a-p , counterpropagating 



lasers with space dependence 9], laser field with heavy nuclei [1014121] and the combination 
of rapidly and slowly varying fields {13I. [l4|. The effect of the temporal laser pulse shape has 
also been investigated [15l4l7]. 



Although the machinery for computin g th e number of pairs produced in the strong exter- 
nal field approximation is well-known [l^1l9|. its evaluation is still a challenging task because 
it is related to a solution of the Dirac equation, which is notoriously hard to solve. For this 
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reason, most of the analytical studies have focused on simple systems. For instance, the pair 
production from a time- varying homogeneous electric field was treated in 5ru| , using differ- 
ent schemes of approximation. The main goal of this article is to revisit this problem from 
a slightly different perspective: the production of pairs from an electric field representing 
counterpropagating lasers in the dipole approximation is evaluated by using a solution of 
the Dirac equation and the adiabatic-impulse model, allowing to study interference effects 
in the pair production spectrum. 

It was recently argued by using a semi-classical approximation that these interference 

se's phenomenon and are responsible for the peak-valley structure 



effects are related to Sto 
in the pair spectrum 1 171 . [20|. The main justification of this result is that interference 
occurs between the semi-classical turning points due to the different phases acquired during 
the time evolution. In this work, this phenomenon is investigated further by using the formal 
analogy of the Dirac equation describing our system with the Driven Two-Level System 
(DTLS). It is well-known that in a certain regime (defined later), the time evolution of the 
DTLS occurs adiabatically, except for some specific times where the system goes through an 
avoided crossing between the "dressed" energy levels. In this case, non-adiabatic transitions 
take place between the lower and upper energy levels. As the avoided crossing region is 
passed many times, these transitions interfere and can lead to constructive or destructive 
interference according to the phase accumulated during the transition and the adiabatic 
evolution. This phenomenon is named the Landau-Zener-Stiickelberg Interferometry (LZSI). 
It is described extensively in {21]], and references therein, and is important in molecular 
physics 22M24I| . In this work, it is shown that the LSZI is also relevant for pair production 
in high intensity laser field and that pair production in lasers proceeds by periodic non- 
adiabatic transitions. To make this connection and for simplicity, the 1-D case is considered, 
which corresponds in 3-D to the production of pairs at zero transverse momentum (p± = 0). 
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20, 



25] 



Other 1-D models have been considered in the literature 

This article is separated as follows. In Section [Til the formalism to compute the rate of 
pairs produced is presented. More precisely, it is shown that the average number of pairs 
produced is related to coefficients in the solution of the Dirac equation. The latter is solved in 
Section [TTT1 in the background field of linearly polarized counterpropagating lasers. In Section 
IIVI a theoretical approach to evaluate the wave function approximately in the adiabatic- 
impulse model is presented. The pair production and numerical results are shown in Section 



\V\ along with an interpretation in terms of nonadiabatic transitions. We take advantage of 
the analogy with the DTLS and evaluate the number of pair produced with the adiabatic- 
impulse model, allowing to understand the spectrum in terms of the LZSI. We conclude in 
Section ED Throughout this work, we use the metric g^ u = diag(l, —1, —1, —1). Also, units 
in which h = c = m = 1 (where m is the electron mass) and e = y/a are utilized in most 
numerical calculations. In this case, the unit length is l u = h/ (mc) ~ 3.86159 x 10 -13 m (38.6 
picometer) while the unit time is t u = h/(mc 2 ) ~ 1.2880885 x 10 -21 s (1.288 zeptosecond), 
as compared to atomic units: Z a . u . = 0.052 nm and t a . u . = 2.4 x 10~ 18 s (2.4 attosecond). 
Note that in these units, the Schwinger field is \eE$\ = 1. 



II. PAIR PRODUCTION FROM STRONG CLASSICAL FIELDS 

The mathematical description of pair production requires a Quantum Field Theory (QFT) 
treatment because it involves particle creation and annihilation. The main tool to calculate 
observable quantities in this framework is a perturbation theory in terms of the coupling 
constant (Feynman diagrams), which allows to evaluate approximately the value of field 
correlators. The latter can be linked to physical observables by using reduction formula 
based on the Lehmann-Symanzik-Zimmermann (LSZ) asymptotic conditions for the field at 



t = =Foo. In these limits, the quantized field operator ^ is known and is given by 26 ] 



*m,out(aO = lim *0M), (2) 

t— >=foo 



dp 



»•< >'»(/') t S -iE^t+ipx , 6 in,out(P) ( s MS 

in,out "m,outU^e -|- . . U mfiut yp)C 



2Ep om ' 2E_ p 



(3) 



(2tt 

where Oin )0 ut(p)) &in,out(p) are annihilation operators that annihilate the "in, out" vacuum as 

ain,out(p)|0in,out) = &in,out (p) \ Oin.out) = °- ( 4 ) 

In these last equations, we have also introduced the asymptotic energies E™ ,out and the free 
positive/negative energy spinors Ui n ,outj ^in,out- The explicit expression of these quantities 
depends on the form of the electromagnetic potential and the gauge chosen; they will be 
described precisely below. On the other hand, the final result for the number of pairs 
produced is gauge invariant. 

When a strong external field such as a laser is involved and coupled to the fermionic 
degrees of freedom, the "naive" perturbation series is no longer asymptotically convergent 



because field insertions, being parametrically of order 1/e, contribute to leading order. The 
field insertions have to be resummed to obtain the leading order contribution. For some 
observables, such as the average number of pair produced from the vacuum (n), this re- 
summation can be performed in the Schwinger-Keldysh formalism by using the Lippmann- 
Schwinger equation 27H29l|. The main result of this procedure is a relation between the 



physical quantity (n) and a solution of the "classical" Dirac equation in the laser back- 
ground field. 

This relation is the starting point of this work and is given by 



27 



29] 



(n) 



dPoutd^in 



lim e iE ™ H ui nt {p)^ q {t,p) 



t—>oo 



(5) 



where ^ q (t,p) is the Fourier transform (with respect to space) of *& q (t, x), the retarded 
solution of the 1-D Dirac equation in coordinate space. Here, the subscript refers to the 
momentum of the initial state: the wave function is subjected to the initial condition 



lim V q (t,p) = v m (q)e iE -^(27r)5(p + q). 

t^t — OQ 



(6) 



27- 



It should be noted here that this condition is derived in the resummation procedure 
291 ] . Physically, this means that the average number of pair produced is computed by 
preparing the system in a negative energy state at t — > — oo, by evolving the wave function 
in time and projecting it onto a positive energy state at t — > oo. Although the Dirac 
equation is solved with an initial state representing a positron, the quantity (n) represents 
the number of pairs produced from the vacuum, without positrons in the initial state. This 
may seem counter-intuitive but it is related to causality and the fact that (n) requires 
the evaluation of Wightman propagators. This can also be understood in the Feynman- 
Stiickelberg interpretation of the positron as an electron evolving backward in time. 

In 1-D, the wave function \1/ is a bi-spinor and the Dirac equation obeyed by ^ q (t,p) is 
then 

id t V(t,p) = [a(cp + A x (t)) + /3mc 2 } V(t,p), (7) 

where m is the electron mass. The Dirac equation is expressed in a gauge where the 
scalar potential Aq = while the vector potential A x (t) is both time-dependent and space- 
independent. Throughout this work, we work in a representation where the Dirac matrices 
are given by Pauli matrices such as a = a z and (3 = <r x . In such representation, Eq. <£T§ is 
supersymmetric [3^] for which semiclassical WKB approximation are exact 31]. 



The covariant measure dp\ 



dp 



is defined with respect to the "asymptotic 



in, out 



(2tt)2E_ 



iin,out 

■± P 



energies 




y/(cp + G in > out ) 2 + m 2 c 4 



(8) 



which are obtained by the on-shell conditions at t = ±00. The G m,ont are constants related 
to the gauge potential as lim^ioo A x (t) = G in ' out . Thus, although the physical electric field 
vanishes asymptotically, the potential may have a non-zero value which depends on the 
gauge chosen. 

The spinors Wi n ,out(p) and tVout(p) are the positive and negative energy solutions of Eq. 
([Z|) with G(t) = G m ' ont and normalized such that ul nout (q)u intOUt (q) = vl nout (q)v- m}OUt (q) = 
2^m,out E X pii c itiy ; they are given by 



As usual, the index in, out represents the limits at t — > =Foo. 

So far, we have defined a general expression for (n). We now specializes this to the case 
of an external laser field. More precisely, the electric field considered is given by 



where E is the field strength and oj is the laser frequency. It represents the field from 
linearly polarized counterpropagating lasers where the space variations are neglected. In 
other words, we consider pair production in the neighborhood of the standing wave anti- 
nodes, where the electric field reaches its maximum value. This field is applied during a 
time interval t G [0, T] where T is the final time. By working in a gauge where the scalar 
potential is zero (Aq = 0), the vector potential is given by 





E(t) = Esm(cot) 



dA x (t) 
dt 



(10) 




where G{t) = 
as F = |e|oE). 



fE^dt' 
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The solution of the Dirac equation in the potential considered can then be written as 



vUq)e lE -^(27i)6(p + q), t e (-oo,0] 

1>(t,p), te[o,T] 

Au ont (p)e- iE ? uH + Bv ont (-p)e iE ™\ t G [T, oo) 



(12) 



where ip(t,p) is a solution of the Dirac equation with the laser vector potential, A, B are 
integration constants that need to be determined from the data at t — T (thus, their value 
depend on T and p) and which allows to have a linear combination of negative and positive 
energy free solutions. Thus, initially, the wave function is given by v- m which represents a 
positron. Substituting the last equation in Eq. (jSj), taking the limit, using the fact that 
u*(p)v(— p) = and vJ(p)u(p) = 2E p , and integrating on the positron momentum (using the 
delta function of the initial state, related to translation invariance), we get 



V f E! ovlt 
HT)) = - dp^-\A(T,p)\\ 



(13) 



where V = 5(0) is the infinite volume. As usual, this diverging quantity is treated by 
redefining (n) as the number of pair produced per unit volume. This is the convention used 
in the rest of this work. 

By requiring the continuity of the solution at t = T, we get the following conditions: 



JE° ut T 



iP(T,p)=A(T,p)u out (p)e-^* i +B (t,p)v ont (-p)e^ 



(14) 



This can be used to compute the constant A which is directly related to pair production via 
Eq. ([TBI . It is a straightforward calculation to obtain 



A(T,p) 



?w,i(p)^i(r,p) + u out)2 {p)^2{T,p) 

2£° ut 



AE° ut T 



(15) 



Thus, we have all the ingredients to calculate the pair production from the vacuum in a 
counterpropagating laser field. In the following, the coefficient A will be evaluated numer- 
ically to obtain the average number of pair produced. The calculation starts by obtaining 
an exact solution of the Dirac equation with the time-dependent background field: the wave 
function ip(T,p) has to be determined. 



III. SOLUTION OF THE DIRAC EQUATION 



Substituting the potential denned in Eq. (fTTj) for t G [0, T] in Eq. ([7]) yields the following 
Dirac equation: 



&z ( C P H cos(cut) ] + a x mc 



(16) 



It should be noted here that this is identical to the equation describing the strongly pe- 
riodically driven two-level system 2jJ, |32] and is also in the supersymmetric form of the 
Dirac equation 30|, ISJj . The formal analogy between the two systems is recovered by letting 
mc 2 — > —A/2, cp — > eo and — — >• A (in the notation of |2l|). This means that for each 
momentum p corresponds a different two-level system. 

Using the explicit expression for Dirac/Pauli matrices, the last equation can be written 
componentwise as 



i-^ I cp H cos(u;£) 

at \ to 



Tpi )2 (t,p) = mc 2 ij 2 ,i(t,p). 



(17) 



These two equations can be decoupled easily to get the following system of differential 
equations: 



dt 2 



UJ 



— iFsm(ut) + I cp H cos(a;t) ) + m 2 c 4 



mc* 



id t — ( cp H cos(o;t) 



jj 



ipi{t,p). 



(18) 
(19) 



Eq. ( TT8l) is a second order Hill's differential equation. It is to be noted that Eq. ( |T8l) contains 
an explicit imaginary part as in optical potential problems, corresponding to absorption 
out of the wave function ipi into the state ip2, i-e. a non-adiabatic transition 33j. As 



demonstrated in 



34j . the last equation can be solved analytically in terms of Heun's function 



(similar equations were also treated in 35 



36]), allowing to evaluate the first component of 



the wave function %p\. The second component can then be found by substituting ip\ into Eq. 
( JT9l) . This way of calculating the solution insures that the general solutions of Eqs. (TT8l) and 
( TT9|) are also solutions of the first order system of equations given in Eq. (|T7|) . at all times. 



To solve Eq. f|T8|) . we follow a strategy similar to [34]: the time domain M + is separated 
into subdomains of length 7t/2uj, parametrized by a positive integer n (see Fig. [1]). Then, Eq. 



( fT8l) is solved for t e At n ■■ 
by matching solutions at points t n 



(2n-l)?r (2n+l)7r 
4ui ' 4u 



the general solution at all times is determined 
^x^j for each interval, using continuity conditions. 



n= (0) 



(1) 



(2) 



(3) 



(4) 



II 1 1 
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V 

in 


1 




1 

7T 

4co 


1 

3tt 
4o> 


1 

4 co 


1 

7n 
4co 

► U 


1 

9tt 
4to 



FIG. 1: The domain M + for t is separated into subdomains parametrized by n G Z + such 



that t G 



(2n-l)w (2n+l)7r 



At the final time t = T, the wave function is matched to the 
free solution u out ,v out . 



The rationale behind this procedure is related to the convergence radius of Heun's functions, 
as will be clarified later. 

The following change of variable is then applied to Eq. f lT8|) for t G At n : 



'7 ij 



(20) 



This transformation maps the time interval At n onto a new variable iz G [— tan(7r/8), tan(7r/8)]. 
The inverse transformation is given by 



2 titt 
— arctanuz) H , 



(21) 



where the arctan function should be evaluated on its principal value. Using this prescription 
and the transformation in Eq. (|20|) on the domain At n , the change of variable is bijective 
and well-defined. A differential equation in terms of the variable z and the model parameters 
is obtained (not shown here for simplicity). The latter is solved by seeking solutions of the 
form: 



exp 



,™2F z 



oo 2 z 2 — 1 



(22) 



where H^\z) is the first linearly independent solution (the second linearly independent 
solution will be denoted by H^jz)). These two transformations convert Eq. (|T8|) into a 



double confluent Heun's equation 



37|: 



d 2 -2z 5 + 4z 3 + a^z 4 -Iz- d z 2 (5^ + (7W + 2a^)z + 5^ 



dz 



(^-1) 3 (^ + 1) 3 



dz 



+ 1) 3 



HV>{z) = 0, 
(23) 



where a^ n \ @( n \ 7^, 8^ n ' are parameters which value depends on the interval considered. 
They are given explicitly by 



a 



7 (») 



8Fcpcos (?f) u - 4iFsin (?f) oj 2 - Aco 2 (c 2 p 2 + m 2 c 4 ) + 2F 2 [(-l) n - 1] 



-16zFcpsin (?f) + 8Fcos (?f) 



CO' 1 



8Fcpcos (?f) 00 - 4iFsin (^) uj 2 + 4w 2 (cV + m 2 c 4 ) - 2F 2 [(-l) n - 1] 



CO' 1 



(24) 
(25) 
(26) 
(27) 



Eq. ( 123|) has singularities at z = ±1 while the point z = is regular. This allows to 
obtain a power series solution around z = with a radius of convergence defined by the 
condition \z\ < 1. Note here that this condition is always fulfilled for all intervals At n 
because iz G [— tan(7r/8), tan(7r/8)| ~ [—0.414,0.414]. This power series is well-known and 
the solution is given by [37J 

Hi n \z) = H D (a^,^ n \^ n \5^,z), (28) 

where Hp is the doubly confluent Heun's function. The partition of the domain has been 
chosen such that the argument of the Heun function z is always within the radius of con- 
vergence, guaranteeing that it is a valid solution and facilitating the numerical evaluation. 
This yields the first particular solution. A second linearly independent solution can be 



found using a well-known procedure |34| : 



Hi n \z) 



— <y(«) * 



'H D (-a (n \(3 in \-f {n \6 (n \z). (29) 

The final result obtained from this is that the first component of the time-dependent wave 
function is given by 



D 



-cr 











(ut 






T, 



(30) 



where A^ n \B^ are integration constants that need to be fixed by initial conditions. The 
second component ^ * s then given by Eq. ( 1T91 . 

The wave function obtained from Eq. ( 1301) is plotted in Fig. [2] for one laser cycle and is 
compared to an accurate numerical solution (the numerical method is described in Appendix 
EJ: both give the same result up to numerical errors. The constants A^ n \ B^ are determined 
by using the continuity of the wave function and its first derivative at points t n . 
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FIG. 2: Real part of the first component of the wave function as a function of time, over 
one cycle of the electric field. The analytical solution is compared to the numerical 
solution. The frequency is u = 0.1, the momentum is p = and the field strength is 

F = 1.0. 



IV. ADIABATIC-IMPULSE MODEL 



Before computing numerical results obtained from the exact solution, it is interesting 
to look more closely at the analogy with the two-level quantum system to gain a better 
understanding of the pair creation process and interference effects. In the adiabatic regime, 
where the laser frequency or photon energy fuv is much less than the gap 2mc 2 and/or the 



magnitude of the vector potential, i.e. when 



38] 



4m 2 c 4 + % > h 2 u 2 , (31) 



CO 2 



it is possible to obtain an accurate approximation of the wave function using the well- 
known adiabatic- impulse model described in [21]. The main advantage of this approach 
is that it yields very simple formula for the transition probabilities such as Landau- Zener 
formula, that allow to obtain new insights into interference phenomena. This model is also 
relevant in our case because the conditions in Eq. (l3Tj) are fulfilled in most prospected 
laser infrastructures aiming at electron-positron production. In a typical laser used to probe 

11 



relativistic effects and pair production, the intensity would be above / ~ 10 24 W/cm 2 , 
leading to an approximate field strength of -Fi aser ~ 2.7 x 10 15 V/m. Using these conservative 
values, the frequency should obey wi aser < 3.5 x 10 18 Hz to fulfill the conditions in Eq. (I3"TI) . 
This implies that in high intensity lasers (with / > 10 24 W/cm 2 ), the production of electron- 
positron pairs can be approximated accurately by the adiabatic-impulse model unless the 
laser frequency is in the 7-ray frequencies. 

Our description of the adiabatic-impulse model starts by considering the adiabatic energy 
of the system given by 



E^(t) = ±J (cp + ^ cos(wt)) + m 2 c 4 . (32) 

The adiabatic energies are plotted in Fig. [3J In this approach, it is assumed that the 
wave function evolves adiabatically at all times except at points where the energy difference 
5-Eadia. (£) = (*) — -^adia (*) * s minimal; at these points, the system undergoes non- 

adiabatic transitions [2l| and the energy difference is 5-E a dia.(^i,2) = 2mc 2 . These points 
correspond to the position of avoided crossings and are given by ut\ = arccos \—*p) and 
ut 2 = 2n — uti (see Fig. [3]). It should be noted here that these conditions are realized only 
if \p\ < F/cu, implying that non-adiabatic transitions do not occur for large momenta since 
then, the adiabatic levels (1321) are well separated at all times. 

In the adiabatic-impulse model, the time evolution is split in two parts (with e a small 
positive time): 

1. Adiabatic evolution: for t G [0, t\ — e], t G [t 1 + e, t 2 — e] or t G [t 2 + e, T]. 

Then, the wave function is given by Eq. flB14[) which is the adiabatic wave function. 

2. Non-adiabatic evolution: for t G [ti — e, t\ + e] and t E [t 2 — e, t 2 + e]: 

Then, the wave Dirac equation can be linearized close to t 12 . The resulting equation 
has a solution in terms of parabolic cylinder function (see Appendix [Bj 

It is then possible to match the wave functions in the two regimes by looking at the asymp- 
totic behavior of each solution at a time t a (which obey t^z *C t a |ii — 1 2 \) and solving for 
the integration constants (see Appendix [B] and [21] for calculation details, and ref. 23] in 
the molecular physics context). The result of this procedure can be casted in a very compact 
notation using transfer matrices. The final result is that (for T G [t 2 , t\ + 2-k/u] ): 

B(T) = U(T, t 2 )NU(t 2 ,t 1 )NU(t 1 ,0)B(0) (33) 
12 



where the vector 



B(t) 



B+(t) 
B-(t) 



(34) 



contains the integration constants of the adiabatic solution. Note that in our case, given the 



initial condition of the wave function in Eq. QI2J, we have that B(0) = (0, ^/2E(0)) T . The 
adiabatic time evolution is generated by the operator 



U(t f ,ti) = exp 



-ia z 



(35) 



while the non-adiabatic evolution is characterized by the matrix iV defined in Eq. ( 1B17I) . 
Using this result, it is straightforward to compute the rate d{n)/dp. We obtain: 



For T G [0,ti]: 



d(n) 
dp 



0. 



(36) 



Initially, the system only has negative energy states and no transition to the positive 
energy states occurs adiabatically. 



For Te [t u t 2 ]: 



d{n) _ j_e; 

dp 



2tt E, 



-P, 



(37) 



When the time reaches t = ti, there is an non-adiabatic transition from the negative 
to the positive energy states with a probability Ps- 



For T <E [t 2 ,ti + 2n/u]: 



d(n) 1 E° v 



out 



dp 2tt E'l 



-AP S {1-P S ) cos 2 ( X + 



(38) 



where \ = J* t * 2 E + (t)dt and is Stoke's phase defined in Eq. (1B18I) . When the time 
reaches t = t 2 , there is another non-adiabatic transition. The wave function coming 
from the negative energy states interfere with the part of the wave function already 
present in the positive energy state, creating an interference pattern characterized by 
cos 2 (x + 4>)- This is the essence of LZSI and is depicted in Fig. |3j 
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FIG. 3: Adiabatic energies in the driven two-level model. The nonadiabatic transitions 
occurs at times ti j2 . The probability of transition is given by Ps while the probability of 
staying in the same state is 1 — Pg. In red are the different paths for transitions from 
negative to positive energy states. After one transition (at t 2 and all times afterwards), the 
part of the wave function in the negative energy states that transits upward, interferes 
with the part of the wave function in the positive energy states. This is the LZSI. 



These features will be seen explicitly in the next section where the pair production will be 
evaluated numerically. 

Note also that these formula can be used to evaluate (n) at later times by applying the 
matrices U and N. For instance, for j laser cycles when T G [t 2 + 2jn /cu, ti + 2(j + l)n/u], 
this would be given by 

B(T) = U(T, t 2 )N [U(t 2 , tjNUitutJN] 3 U(t 2 , h)NU(ti, 0)B(0), (39) 



for j G Z + . An explicit expression of the matrix [U(t 2 , ti)NU(ti, t 2 )N] J is given in 21]. 



V. NUMERICAL RESULTS 

In this section, the rate of electron-positron pair production is calculated numerically 
using the exact solution and the adiabatic-impulse model. The first result concerns the 
quantity d(n)/dp in the adiabatic regime, which is plotted in Fig. HI In this figure, we 
also include the position where non-adiabatic transitions take place. It is clear from this 
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d(n)/dp 




FIG. 4: Rate of pairs produced from the counterpropagating laser. The frequency is 
uj = 0.1 and the field strength is F = 1.0, which insures that the system is in the adiabatic 
regime. The black line shows the position in (p, T)-space where the non-adiabatic 

transitions take place. 

figure that qualitative changes occur at these points. This can be understood very clearly by 
looking at the theoretical results obtained from the adiabatic-impulse model, in Eqs. ( 1361) to 
(138p . In the first instants, there is no pair production because the system starts in a negative 
energy states and there is no transition to the positive energy states when the wave function 
evolves adiabatically. For a given momentum p, when the time reaches t = t\, there is a 
non-adiabatic transition and pairs start to be produced with a rate given approximately by 
Eq. ( |37l) . Later in the time evolution at t = £2, a second transition happens and interferes 
with the preceding one, resulting in an interference pattern described by Eq. (1381) . This is 
the well-known LZSI. After this, each time the system crosses a non-adiabatic transition, a 
part of the negative energy states traverses to the positive ones and a different interference 
pattern emerges. The corresponding average number of pairs produced (the spectrum is 
integrated on p at each time) in shown in Fig. [5j 

These results show a very good qualitative agreement between the two theoretical ap- 
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FIG. 5: Average number of pairs produced and electric field E(t). The frequency is 
oj — 0.1 and the field strength is F = 1.0. 



proaches. To investigate this comparison more quantitatively, the pair production spectrum 
after one laser cycle (for T = 2tt/uj) is plotted in Fig. |6l along with the prediction of the 
adiabatic- impulse model (more precisely, the envelope obtained from Eq. (138]) ). The spec- 
trum shows the characteristic peak- valley structure of an interference pattern, in agreement 



with the results obtained in [16 
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20(. The latter is well described by Eq. ( 1381) (it was 



also verified that the maxima and minima of the spectrum corresponds to those of Eq. (I38j) ) 
and thus, this effect is due to the LZSI. Also, it should be noted that again, both theoretical 
approaches yields very similar results. 

Finally, in Fig. [TJ the rate of pairs produced is presented for a smaller value of field 
strength (F = 0.1). In this case the magnitude of the interference patterns are decreased 
significantly and the spectrum is peaked at the value of the non-adiabatic crossings. The 
most likely explanation for this behavior is related to the non-adiabatic transition time, 



which can be estimated as 



21 



38 



39| 



t n .-a. ~ 



mc 



(40) 



The adiabatic-impulse model requires that t n .- a . "C vr/w (the transition time should be 
much shorter than a half-cycle) to make sure that each transition is independent and well 
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FIG. 6: Spectrum of pairs produced at T = 2n/cu. The frequency is uj — 0.1 and the field 
strength is F = 1.0. The exact result is compared to the envelope obtained from the 

adiabatic-impulse approximation. 

separated from each other in time. Clearly, for F — u — 0.1, this condition is not fullfilled 
and the impulse-model is not valid for these parameter values. 



VI. CONCLUSION 



In this article, the production rate of electron-positron pairs from high intensity linearly 
polarized counterpropagating lasers has been considered. An exact solution of the Dirac 
equation in terms of Heun's function and numerical methods has been used to compute this 
observable. The results has been compared to the ones obtained from another theoretical 
approach called the adiabatic-impulse model. The latter is based on the adiabatic approxi- 
mation and allows to obtain simple expressions for the wave function. The results obtained 
from both methods were very consistent with each other in the adiabatic regime and when 
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d(n)/dp 




FIG. 7: Rate of pairs produced from the counterpropagating laser. The frequency is 
u = 0.1 and the field strength is F = 0.1. The qualitative behavior is different from the 
adiabatic regime. The black line shows the position in (p, T)-space where the 
non-adiabatic transitions take place. 



the non-adiabatic transition time is much shorter than a half-cycle. Therefore, it has been 
concluded that the adiabatic-impulse model is an accurate theoretical tool that may be used 
for other more complex systems. It was demonstrated that pair production occurs through 
periodic non-adiabatic transitions and that these transitions resulted in a complex interfer- 



ence pattern in the pair spectrum. This is very similar to the results obtained in [16|, [17|, 20 ]. 
This phenomenon has been related to the well-known LZSI by using the formal analogy of 
our Dirac equation with the DTLS. 



Appendix A: Numerical method 

The numerical method used in this work is inspired from the spectral methods developed 



m 



40 



431 ]. where an operator splitting scheme is used in momentum space to evolve the 



solution in time. Other approaches where the Dirac equation is solved in momentum space 
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r ( F 




a z [cp+ — 







can be found in 44j, |45( . The solution of the Dirac equation in Eq. (Tl6]) can be written 
formally as 

^*(0,p) (Al) 

where "T" stands for the time-ordering operator. This expression can be used to solve the 
Dirac equation numerically in momentum space. First, the time domain is separated into 
TV time increment having a size St. Then, it can be shown that the last expression can be 
approximated by 



tt(t, p) = W(t, t N )W(t N , tpf—x) ■ ■ ■ W{h, t )*(0, p) + 0((5t) 3 ) 



(A2) 



where L- = jSt and the evolution operators are given by 



(a z cp + a x mc 2 )5t + a z I dt! — cos(out') 

Jti-! ^ 



(A3) 



The last equation can be computed explicitly by using the properties of Pauli matrices. It 
can then be written as: 



ITT/, , \ TT / \ -O'X^X ~\~ Q-Z&Z ■ / \ 

W(ti, ti-i) = h cos(a) — i sm(a) 



(A4) 



where 



a x = mc 2 5t, 

F 

a z = cpSt H — - sin(a;ti) - sin(a;tj_i) , 
a = \Ja 2 £ + a 2 z . 



This results in a numerical method for which the error is 0(St 3 ) 46] 



(A5) 
(A6) 
(A7) 



Appendix B: Solution at the non-adiabatic transition 



In this appendix, the transfer matrix around a non adiabatic transition is derived, fol- 



lowing the work exposed in 



2J24 



24) (and references therein). 
In the neighborhood of t± 2, the vector potential can be linearized in if and the resulting 
Dirac equation is given by 



id t 'ip{t',p) = \^co z vt! + a x mc 2 ] ip{t',p), 



(Bl) 
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for tx t 2, respectively, and where i! = t—ti >2 (in the following, we suppress the prime notation). 

/ 2 2 2 

Here, we have v = Fyl — c ^ 2 P . The last equation is formally equivalent to the Landau- 
Zener transition which has a well-known solution in terms of parabolic cylinder function. 
The latter can be found by writing the last equation componentwise and by decoupling the 
two resulting equations. This yields 



d 2 

u i 2.2 -i 2 4 

— + v t — iv + m c 
at 2 



d 

l dt +Vt 



The solution of Eq. (jB2[) is found by a change of variable given by z 
transforms the equation to 



r ai 



dz 2 ~I~2~ l 



0. 



(B2) 
(B3) 

2ve l ^t, which 
(B4) 



where we defined 5 = ^y 2 -- The last equation has a solution given by [4 



(B5) 

where C\ y2 are integration constants and D u (z) is the Whittaker parabolic cylinder function. 
It is then a straightforward calculation to obtain the second component using Eq. (IB3P and 
the recurrence relation of the parabolic cylinder functions. We get that 



2ve~ l -t 



(B6) 



lim^(*)~C 2 -^-e-^), 
t^t a -1(1 + 



The next step is the evaluation of the wave function far from the non-adiabatic transition 
region, that is when t = \t a \ ^> 1. Using the asymptotic expansions for D u , we obtain 

(B7) 
(B8) 
(B9) 
(BIO) 



lim V> 2 (*) ~ [-Cie* 6 + C 2 e-i 6 ] 



?<5 , n e 4 4 r -i<Ht) 



t-¥t 



V6 



lim Mt) ~ C 

t-f-ta 



fSA<f,(t) 



-e 4 e 



1 r(i + 2<5) 

lim fo(t) ~ [-de^ 5 + C 2 e^ 5 ] 



e 4 4° 



-i<t>{t) 



V6 



where we defined the time-dependent phase as 

vt 2 



0(t) 



5 In (V2vtj . 



(Bll) 
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We would like to match these asymptotic wave functions to the adiabatic wave function far 
from the transition times tip. The adiabatic wave function is obtained as follow. 

In the adiabatic approximation, the wave function looks like 

-t 



^adia.(*) = f eX P 



=Fi / E(t')dt' 



(B12) 



where the ± denotes positive and negative energy solution, respectively, and ^ are adiabatic 
coefficients to be determined which obey Ic^y? 1 * 1 ! <C \E{t)(p ± \ for all times. By substituting 
this into the Dirac equation 



idtif>(t) = [a z P(t) + a x mc 2 ] i/>(t), 



(B13) 



where P(t) is the canonical momentum for the system under consideration, normalizing the 
wave function such that iV^ia | 2 = 1, we arrive at the following general solution (which is a 
linear combination of positive and negative energy solutions): 

-t 



V'adia.(t) 



± 



B ip exp 



Ti / E(t')dt' 



7T 



4 



(B14) 



where 



/E+P(t) 




/E-P(t) 


V 2E 


; <p = 


V 2E 


l /E-P(t) 


/e+p® 


LV 2£ J 




L V 2E J 



and where are integration constants. Note that the factor n 



(B15) 



74 in the phase appears 



when the next order in the adiabatic approximation is considered 221 ] . 

Here, we are considering two times ±t a which are lying in a region close to the transition 
region but which are much longer than the typical non adiabatic transition time tj^z, that 
is £lz *C \t a \ \ti — ^2 1 - In this region, the time evolution of the wave function can be 
described accurately by the linearized Dirac equation, Eq. fIBlj) (thus, Pit) = —vt), while 
still being far from the transition times. In this case, assuming that v\t\ 3> mc 2 also holds, 
the adiabatic solution can be simplified to give 

BTe ^(t a )-i|[ln(5)-l]-i| 

±_B± e -^(tu)-HfM*)-i]+*3 



^adia.(±^a) 



(B16) 



The time dependence of the last expression is the same as the asymptotic solutions in Eqs. 
( IB7j) to (IBlOj) . allowing to match the solutions at t = =)=t a and thus, to determine the transfer 
matrix that allows to link the solution at negative times to the one at positive time. Solving 
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for the integration constants in Eq. (1B16j) . we find that the non-adiabatic transition can 
be characterized by the following time-independent transfer matrix 2l|, |24j : 



N 



VP's 



-VP's 



(B17) 



where the Stoke's phase is defined as 



7T 



- + 5[\n(5) - l]+argr(l-i<y). 



(B18) 



The quantity P$ = e~ 2nS is Schwinger's result for the pair probability creation in a constant 
field (note that it is also the Landau-Zener transition probability). The matrix N connects 
the wave function before and after the non-adiabatic transition. 
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